function sdm

%  SDM for 2x2 system
%        [a b][x1] = [b1]
%        [c d][x2] = [b2]

% clear all previous variables and plots
clear *
clf

% pick an example (iexample=1,2)
iexample=2;
tol=0.1;

% set parameters
if iexample==1
	a=2; b=2; c=b; d=3; b1=1; b2=-1; 
%	xL=0; xR=6; yL=-5; yR=1;
	xL=0; xR=4; yL=-6; yR=2;
	cxL=3; cxR=10; cy=-2;
elseif iexample==2
	a=5; b=4.99; c=b; d=5; b1=1; b2=-1; 
	xL=0; xR=80; yL=-80; yR=10;
	cxL=0; cxR=100; cy=-100;
end;
x=[1, 1];

% generate surface data for plotting
xd=linspace(xL,xR,60);
yd=linspace(yL,yR,60);
[X,Y]=meshgrid(xd,yd);
Z = 0.5*(a*X.^2 + 2*b*X.*Y + d*Y.*Y) - b1*X - b2*Y;
%surf(X,Y,Z)

% determine elevations to draw contours
cx=linspace(cxL,cxR,18);
for i=1:18
	cv(i)=0.5*(a*cx(i)^2 + 2*b*cx(i)*cy + d*cy*cy) - b1*cx(i) - b2*cy;
end;

% get(gcf)
set(gcf,'Position', [1034 527 571 493]);
plotsize(571,493)


% draw countors and put in labels
contour(X,Y,Z,cv,'b')
hold on
box on
xlabel('v_1-axis','FontSize',14,'FontWeight','bold')
ylabel('v_2-axis','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',14);

% have MATLAB use certain plot options (all are optional) 
if iexample==1
	axis([0 4 -3 1])
	set(gca,'ytick',[-3 -2 -1 0 1]);
	axis square
elseif iexample==2
	axis([0 80 -80 1])
	set(gca,'ytick',[-80 -60 -40 -20 0]);
	axis  square
end;

% start iteration
xpoints=x(1);
ypoints=x(2);
r=[b1-a*x(1)-b*x(2),b2-c*x(1)-d*x(2)];
diff=[1,1];
counter=1;
error=norm(diff);
while norm(diff)>tol & (counter<30)
	counter=counter+1;
	q=[a*r(1)+b*r(2),c*r(1)+d*r(2)];
	alpha=dot(r,r)/dot(r,q);
	diff=alpha*r;
	x=x+diff;
	r=r-alpha*q;
	xpoints=[xpoints, x(1)];
	ypoints=[ypoints, x(2)];
	plot(xpoints,ypoints,'r-','linewidth',1)
	plot(xpoints,ypoints,'o','linewidth',1.1,'markeredgecolor','r','markersize',8)
	error=norm(diff);
	title(['iteration number = ',num2str(counter),'    Iteration Error = ',num2str(error)],'FontSize',14,'FontWeight','bold')
	pause
end;

% plot(xpoints,ypoints,'r-','LineWidth',1)
% plot(xpoints,ypoints,'o','LineWidth',1.1,'MarkerEdgeColor','r','MarkerSize',8)

hold off

%A=[[a b];[c d]];
%bb=[b1;b2];
%exact=A\bb
%error=norm(exact-x');
%fprintf('\n  Iterations= %i  Error =  %e  %e \n\n',counter,error);
%eigenvalues=eig(A)


% subfunction plotsize    '
function plotsize(width,height)
siz=get(0,'ScreenSize');
bottom=max(siz(4)-height-95,1);
set(gcf,'Position', [2 bottom width height]);